CN106646645B - 一种重力正演加速方法 - Google Patents

一种重力正演加速方法 Download PDF

Info

Publication number
CN106646645B
CN106646645B CN201611252162.XA CN201611252162A CN106646645B CN 106646645 B CN106646645 B CN 106646645B CN 201611252162 A CN201611252162 A CN 201611252162A CN 106646645 B CN106646645 B CN 106646645B
Authority
CN
China
Prior art keywords
gravity
point
cell
layer
octree
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
CN201611252162.XA
Other languages
English (en)
Other versions
CN106646645A (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 CN201611252162.XA priority Critical patent/CN106646645B/zh
Publication of CN106646645A publication Critical patent/CN106646645A/zh
Application granted granted Critical
Publication of CN106646645B publication Critical patent/CN106646645B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • Remote Sensing (AREA)
  • Computer Graphics (AREA)
  • Theoretical Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Medicines Containing Antibodies Or Antigens For Use As Internal Diagnostic Agents (AREA)
  • Pharmaceuticals Containing Other Organic And Inorganic Compounds (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种重力正演加速方法，包括以下步骤：步骤1、将三维地质体剖分成多个四面体单元；分别对观测点和三维地质体建立对应的八叉树结构；步骤2、将三维地质体划分为近区与远区，步骤3、分别计算近区与远区重力参数；基于近区所包含的四面体单元的重力积分无奇异性解析表达式直接计算近区的积分，得到近区重力参数；并基于观测点八叉树和源八叉树，采用自适应快速多极展开算法计算远区的积分，得到远区重力参数；步骤4、将步骤3中计算得到的近区重力参数与远区重力参数求和，得到观测点重力参数。本发明能够基于非结构化网格剖分，对任意分布的观测点进行重力正演计算，大幅度地提高了重力正演计算的效率。

Description

一种重力正演加速方法
技术领域
本发明涉及一种地球物理领域的重力正演方法，特别涉及大规模航空重力、卫星重力等的高精度、快速正演方法。
背景技术
重力勘探是地球物理勘探的一种重要的方法，并被广泛应用到矿产资源勘探、地壳结构与莫霍界面探测和环境地球物理勘查中。测量的重力数据与地下密度分布体有着直接的关系，因此，通过对地下密度进行反演可以用来进行地下密度异常体的判别。随着重力勘探技术的发展，航空重力技术的逐渐成熟，对于给定密度分布的目标体，重力勘探往往含有较多的观测点，传统的重力正演方法计算耗时的特征越来越明显。因此，高效的正演显得尤为重要，目前，常用的重力计算加速算法主要包括两种：基于快速傅里叶变换(FastFourier Transformation,FFT)的重力正演和基于小波变换(wavelet transformation)的重力快速正演。以上两种方法计算效率快、精度高，但是都需要对三维地质体进行均匀的六面体网格剖分，这使得网格剖分不能很好的拟合复杂模型，方法不适合处理带地形的三维问题，同时该方法只能计算观测点在同一个平面上的重力正演。非结构化网格能够拟合复杂模型的优势日益突出，被广泛应用在三维数值模拟中。
因此，有必要设计一种高效、高精度的基于非结构化网格剖分的观测点任意分布的重力正演方法。
发明内容
本发明所解决的技术问题是，针对现有技术的不足，提供了一种重力正演加速方法，能够基于非结构化网格剖分，对任意分布的观测点进行重力正演计算，大幅度地提高了重力正演计算的效率。
本发明的技术方案为：
步骤1、模型建立、区域离散与八叉树建立：
如图1所示，三维地质体Ω的密度分布为ρ(r′)，在观测点r处产生的重力位φ(r)、重力梯度g(r)和重力梯度张量Tg(r)可以表示为：
其中，为梯度运算符，R＝|r-r′|为观测点r到源点r′的距离；dv表示体积分，ρ(r′)为三维地质体Ω在源点r′处的密度，G为引力常数，G＝6.674×10-11m3kg-1s-2；
模型建立：根据观测点与三维地质体(又称为源)的特征，设计三维地质体模型参数(模型坐标及密度信息)和观测点参数(观测点坐标)。
区域离散：对三维体地质体模型采用四面体单元进行网格离散，将三维地质体离散成多个四面体单元Ti,i＝1,2,...,N，N为四面体单元的个数。因此，三维地质体Ω可表示为：
八叉树建立：八叉树是一种用来描述三维空间的树状数据结构，本发明采用八叉树的数据构建方法，实现对观测点以及地下地质体建立相应的八叉树数据结构，即观测点八叉树和三维地质体单元(源)八叉树，如图3所示；
对源建立八叉树：构建一个立方体单元包含三维地质体，并将该立方体单元称第0层细胞；然后以第0层细胞为母细胞，将其细分成8个子细胞，该8个子细胞称为第1层细胞；仍然，按照上述剖分方式对每个子细胞进行细化，后续每剖分一层得到的子细胞依次记为第j(j＝2,3,4…,q)细胞层子细胞，直到每个子细胞里包含的四面体地质单元数少于设定的个数为止。最后一层的子细胞称为叶子细胞，同时删除不含有四面体地质单元的细胞。
同样地，构建观测点的八叉树数据结构：构建一个立方体单元包含所有的测点，称为第0层细胞；然后以第0层细胞为母细胞，将其细分成8个子细胞，该8个子细胞称为第1层细胞；仍然，按照上述剖分方式对每个子细胞进行细化，后续每剖分一层得到的细胞依次记为第i(i＝2,3,4…,p)层细胞直至每个细胞包含的观测点数目少于设定的个数为止，同时删除不含有观测点的细胞，最后一层的子细胞称为叶子细胞，。
步骤2、近区与远区的分离：
常规的重力计算方法是每个观测点对每个小的四面体地质单元进行一次积分，然后叠加求和得到该测点的重力值。快速多极展开算法的目的就是通过将积分区域划分成远区和近区，近区积分与常规的重力计算方法一致，对于远区的积分，通过选取扩展点，即每个细胞的中心点，将每个观测点与所有的四面体地质单元的积分转换成观测点所在细胞的中心点与四面体单元所在细胞的中心点的连接，从而减少了积分计算的次数，达到加速的目的。
为了达到这一加速的目的，将观测点与三维地质体之间的连接分为三种连接关系：第一种为三维地质体子细胞中心点与母细胞中心点的连接，通过M2M转换连接；第二种为三维地质体细胞中心点与观测点母细胞中心点的连接，通过M2L转换连接；第三种为观测点子细胞中心点与母细胞的中心点的连接，通过L2L转换连接。如图3所示。
从重力位的计算公式可以看出，积分核函数含有1/R项，当R趋近于0时，积分含有奇异性；当观测点r与三维地质体Ω内部积分点r′的距离趋近于0时，即R趋近于0，因此，将含有奇异性积分的区域称为近区，反之，其他区域称为远区。但是，实际计算过程中，对于每一个观测点r，若三维地质体Ω中某个四面体单元的中心点到观测点r的距离小于某一设定的数值ε【ε趋近于0，可取1-4个四面体单元的长度】，则认为该四面体单元所在区域为近区，记为Tnear；其他区域为远区，记为Tfar。如图2所示。那么，区域Ω可以表示为：
Ω＝Tnear+Tfar (2)
步骤3、分别计算近区与远区重力参数：
将重力正演的积分划分为对近区的积分和对远区的积分；基于近区所包含的四面体单元的重力积分无奇异性解析表达式直接计算近区的积分，得到近区重力参数；并基于观测点八叉树和源八叉树，采用自适应快速多极展开算法计算远区的积分，得到远区重力参数；然后通过求和得到重力正演结果。
近区重力参数包括在观测点r处产生的近区重力位、近区重力梯度和近区重力梯度张量，分别记为和远区重力参数包括在观测点r处产生的远区重力位、远区重力梯度和远区重力梯度张量，分别记为和
步骤4、将步骤3中计算得到的近区重力参数与远区重力参数求和，得到观测点重力参数，包括在观测点r处产生的重力位φ(r)、重力梯度g(r)和重力梯度张量Tg(r)，其计算公式分别为：
以下对所述采用自适应快速多极展开算法计算远区的积分，得到远区重力参数的原理和步骤进行详细说明：
为了使用快速多极展开算法，需要对积分点进行转换，即对积分核函数1/|r-r′|进行展开，当满足条件时，积分核函数可以展开为：
其中，为源点r′所在叶子细胞的中心点(r′的扩展点)，Rnm和Snm分别为常规的(R)和奇异的(S)球谐函数，为Snm的复共轭。这两种球谐函数可以表示为：
其中，(r,θ,φ)为点x的球坐标表示，r>0，0≥θ≥π，0≥φ≥2π。为伴随勒让德函数，定义为：
其中，Pn为标准的n阶勒让德多项式。
此时，将积分核函数分成了两部分，一部分是关于观测点r和的函数另一部分为关于r′和的函数目的是为了将观测点r从积分中分离出来。那么，重力位可以表示为：
其中，为重力位在rc p处的多极距，叶子细胞的多极距直接通过多极距公式计算。多极距只与r′和有关，与观测点r无关，不随观测点的变化而变化，这就意味着，即使有很多个观测点，多极距也只需要计算一次。如果多极距一旦被计算好，当计算观测点的远区积分时，只需要进行一些简单的代数操作，即计算依赖于观测点的基础函数与不依赖于观测点的多极距的乘积。对于大规模重力计算，该方法可以大幅度的减少计算时间。如果所有的源细胞都在远区，那么计算M个观测点的φ(r)far的时间复杂度将从O(NM)减少为O(bM+N)，其中，N为四面体单元的个数，b为细胞数。为了达到进一步减少计算时间的目的，需要借助自适应多层的方法，此时，时间复杂度可以减少为O(MlogN)。
设点rc i为源八叉树第i层细胞的中心点，点为源八叉树第(i+1)层细胞的中心点，点rc i所在细胞为点所在细胞的母细胞；其中i＝0,1,2…,p-1；最高层细胞即第0层细胞的中心点为rc 0；对于母细胞，通过M2M转换公式，对该母细胞的所有子细胞累加求和得到母细胞的中心点处的多极距；
首先，通过点rc i+1处的多极距与以下M2M转换公式计算点rc i处的多极距
然后，通过以下公式计算远区重力位
母细胞的多极距通过其子细胞的多极距和M2M转换公式计算，M2M转换公式实现了子细胞与母细胞之间的联系，同样Mnm(φ,r′,rc i)包含两部分，一部分为与rc i+1和rc i有关的Rn′m′(rc i+1-rc i)，另一部分为与r′和rc i+1有关的即为子细胞的多极距。可以看出，通过M2M转换，点rc i处的多极距Mnm(φ,r′,rc i)只需要进行简单的几何运算即可求得，从而可以快速求得远区重力位。
类似地，如果点为观测点的局部扩展点，当满足条件时，积分核函数可以展开为：
那么，重力位在远区的积分可以表示为：
其中，为已知的基础函数，表示处的局部系数，为未知的局部系数，可以通过计算远区Tfar的积分得到。但是，也可以通过已经计算好的多极距来计算局部系数从而达到减少计算量，加快计算速度的目的。
借助源细胞中心点的多极距来计算扩展点的局部系数具体计算方法为：
首先，对进行展开：
这里，此时，将展开成两部分，一部分为与r′和相关的另一部分为与和相关的可以看出，与观测点及其扩展点无关，很好的将观测点与源点分离开。由此得到M2L转换公式为：
由以上公式求得处的局部系数。
设点为观测点八叉树第(j+1)层细胞的中心点，点为观测点八叉树第j层细胞的中心点，点所在细胞为点所在细胞的母细胞，j＝0,1,2…,q-1；通过点处的局部系数和以下L2L转换公式计算点处的局部系数
那么，点的局部系数可以通过简单的代数求和直接计算得到。使用上述L2L转换，可以最大程度的减少计算时间。
然后，通过以下公式计算远区重力位：
逐层计算，从上往下，便可以计算得到每个观测点的重力位的远区积分。
基于同样的原理，可以快速计算远区重力场和远区重力梯度张量：
通过以下公式计算远区重力场和远区重力梯度张量：
用于计算远区重力场和远区重力梯度张量的局部系数与计算远区重力位的局部系数相同。这就意味着，如果计算远区重力位时，点处局部系数通过L2L，M2L和L2L转换公式计算得到，只需要再计算局部系数和局部基本的函数和的乘积，便可以计算得到重力场和重力梯度张量。
计算局部基础函数和需要使用坐标转换，这里xi(i＝1,2,3)表示三个坐标分量，表示三个方向的单位向量。那么，可以写为：
可以写为：
其中，yi(i＝1,2,3)表示为球坐标系下三个分量，y1＝r>0，y2＝θ和0≤θ≤π，
xi和yi之间的关系可以表示为：
所述步骤3中，采用基于近区所包含的四面体单元的重力积分无奇异性解析表达式直接计算近区重力参数的计算公式为：
其中，R＝1/|r-r′|。
本专利采用基于四面体单元的积分解析计算的方法对上述几种积分进行无奇异性、高精度的解析计算。
对于任意的四面体单元T，借助矢量恒等式，
标量形式可以重新写为：
其中，为四面体单元T的第i个三角面，di＝ni·(r′-r)表示向量(r′-r)在面的法向ni上的投影。当点r靠近点r′时，面积分ki包含对1/R的含有弱奇异性的积分。
借助梯度定理，向量形式可以表示为：
其中，ki定义和前面一致，因此，如果ki计算了，那么体积分A和D便可以计算得到。
可以表示为：
下面，给出面积分和的解析表达式。下面使用K表示三角面用k和k-1分别表示含有奇异性的面积分和这里，|r-r′|表示点r到点r′的距离。
如图5所示，点o为测点r在面K的投影点，并设点o为局部极坐标系的原点。为面K的面外法向，那么，有：
其中，ω0可以为任意实数，ω0＝0表示点r在面K上。为极坐标系下的向量，从点o指向积分点r′，可以写为其中单单位向量，ρ≥0为向量的长度。此时，R＝|r-r′|可以表示为 的奇异性取决于ω0和ρ的数值，当ω0＝0时，在点o附近的无穷小得区域里积分存在奇异性，因此，将积分区域K分为两个部分：可能含有奇异性的无穷小区域C(ε)和其他不含奇异性的区域K-C(ε)，那么有：
K＝[K-C(ε)]+C(ε) (25)
无穷小区域C(ε)为一个角度为α半径ε→0的圆形区域。那么α可能为：
借助如下两个向量恒等式(Wilton,et al,1984)：
其中，为对点r的梯度，表示对积分点r′的面梯度。那么，有：
在积分区域K-C(ε)，R-1无奇异性，借助梯度定理：
其中，为在边上的投影在外法向上的投影，在边上，为常数。类似地，
其中，积分是与点o相关的奇函数。积分只有当观测点在三角面上时才消失，此时a(o)＝2π。当观测点在三角面的边上或者点上时，积分存在奇异性。因此，有：
现在，考虑两个一维线积分。建立一维局部坐标系，定义：
其中，为边的切向方向， 为面K的面外法向，为边的外方向。s＝0表示向量与边垂直，其值可以表示为：
这里，和可以为正数，也可以是复数，表示点o在边的延长线上，表示点o在边上。现在，距离R可以表示为这里等于向量的方向与方向的点积。 和为点r到点和点的距离。和为点r到边的距离。
现在，一维线积分可以表示为：
那么，
这里，角度a(o)可以表示为：
当观测点在边上时，将使数值计算不稳定，为了避免该问题，引入此时，有：
其中，此时，面积分没有奇异性。
从上可知，其中的第二个积分有如下结果：
当时，第一部分可以写为(Gradshteyn,et al,2014)：
当时，如果积分存在奇异性，从新计算积分
将近区与远区对观测点r处产生的重力位、重力梯度和重力梯度张量求和得到整个区域对观测点r处产生的重力位、重力梯度和重力梯度张量，其计算公式分别为：
有益效果：
常规的重力计算方法是每个观测点对每个小的四面体地质单元进行一次积分，然后叠加求和得到该测点的重力值。本发明采用快速多极展开算法，通过将积分区域划分成远区和近区，近区积分与常规的重力计算方法一致，对于远区的积分，将每个观测点与所有的四面体地质单元的积分转换成观测点所在细胞的中心点与四面体单元所在细胞的中心点的连接，从而减少了积分计算的次数，实现加速的目的。为了达到这一加速的目的，分别建立观测点与三维地质体八叉树结构，将观测点与三维地质体之间的连接分为三种连接关系：第一种为三维地质体子细胞中心点与母细胞中心点的连接，通过M2M转换连接；第二种为三维地质体细胞中心点与观测点母细胞中心点的连接，通过M2L转换连接；第三种为观测点子细胞中心点与母细胞的中心点的连接，通过L2L转换连接。计算三维地质体八叉树中各个细胞的多极距或观测点八叉树中各个细胞的局部系数，从而计算远区重力位；三维地质体八叉树中叶子细胞的多极距直接通过多极距公式计算，母细胞的多极距通过M2M转换公式计算；多极距与观测点r无关，不随观测点的变化而变化，这就意味着，即使有很多个观测点，多极距也只需要计算一次。如果多极距一旦被计算好，当计算新的观测点的远区积分时，只需要进行一些简单的代数操作。对于大规模重力计算，该方法可以大幅度的减少计算时间。进一步地，本发明也可以通过计算观测点八叉树中各个细胞的局部系数，从而计算远区重力位；观测点八叉树中最高层细胞的局部系数通过三维地质体八叉树中最高层细胞的多极距和M2L转换公式计算，再通过L2L转换公式计算子细胞的局部系数。
本发明可以进行大规模重力高精度、快速正演计算与大规模重力地形改正计算，大幅度地提高了计算的效率。
附图说明
图1为典型的任意复杂的重力勘探示意图，r表示观测点，r'表示三维地质体Ω内的任意一点，三维地质体的密度分布为ρ(r′)，h表示观测点相对地面的高度。
图2为三维体Ω分成近场区Tnear和远场区Tfar示意图。
图3为观测点八叉树和源八叉树中M2M、M2L和L2L转换示意图。
图4为基于快速多极展开算法的重力场计算方法的流程图。总场通过对近区场和远区场求和得到，近区场直接采用解析公式计算，远区场通过快速多极展开算法计算。
图5为基于四面体单元无奇异性积分局部坐标系示意图；
图6为六面体模型示意图。六面体的密度为2000kg/m3，并被剖分成一系列的四面体单元，观测平面距离六面体的上表面的距离为h。
图7为本发明的基于四面体单元的解析表示式计算结果与基于六面体单元的解析表达式计算结果的对比图。图(a)为基于六面体单元的解析表达式计算结果；图(b)为本发明的基于四面体单元的解析表示式计算结果；图(c)为二者的相对误差图。
图8为本发明的基于快速多极展开算法的重力场计算结果与基于四面体单元的解析表达式计算结果的对比图。图(a)为本发明的基于四面体单元的解析表达式计算结果；图(b)为本发明的基于快速多极展开算法的重力场计算结果；(c)为二者的相对误差图。
具体实施方式
以下结合附图和具体实施方式对本发明作进一步的说明。
本发明涉及的重力计算方法包括以下步骤：
步骤1、测点和源文件的设计：确定测点的个数与坐标、密度体位置坐标文件。并进行非结构化网格剖分。分别对观测点与密度体建立对应的八叉树结构，将积分区域分为近区与远区。
步骤2、远区场计算：通过自适应快速多极展开加速算法计算远区场积分。
步骤3、近区场计算：通过解析表达式计算高精度、无奇异性的计算近区积分。
步骤4、计算观测点重力场：通过对近区场积分和远区场积分求解计算得到观测点重力场值。
以下为本发明计算一个六面体重力场的实例。
如图6所示，设单一六面体的密度为2000kg/m3，长100m，宽100m，高100m。令六面体的上表面中点为笛卡尔坐标系的原点，向下方向为z轴的正方向。M＝210201个观测点分布在六面体上表面上方高度为h＝100m的平面上，观测平面为200m×200m。采用TetGen将六面体剖分成N＝20415个四面体单元。
为了验证本发明方法的正确性，进行一下两个验证：四面体重力场的解析表达式的正确性验证和采用快速多极展开的重力场计算的正确性验证。首先验证四面体重力场的解析表达式的正确性验证，将本发明的解析解计算的结果与六面体单元的重力解析解计算结果(Li,et al,1998)进行对比，只对比了重力位与垂直方向的重力梯度两个结果，如图7所示，两者的结果十分吻合，重力位的最大误差为6.0×10-11％，垂直方向的重力梯度的最大误差为4.0×10-9％。说明本发明推导的四面体重力场解析表达式的正确性。下面验证采用自适应快速多极展开算法的重力场计算的正确性，将结果与采用四面体解析表达式计算的结果作对比，同样只比较重力位与垂直方向的重力梯度两个分量，如图8所示，重力位的最大误差为1.0×10-5％，垂直方向的重力梯度的最大误差为5.0×10-9％，二者的结果十分吻合，说明本发明基于快速多极展开算法的重力场计算方法的正确性。使用基于四面体单元的解析表达式进行计算的方法计算时间为63.6s，使用本发明基于快速多极展开算法的重力场计算方法的计算时间为7.1s，加速比为8.95倍。并且，观测点数越多，加速比越大。同时采用OpenMP并行计算加快计算，当M＝1002001时，加速比达到8661倍。
参考文献
Li X,Chouteau M.Three-dimensional gravity modeling in all space[J].Surveys in Geophysics,1998,19(4):339-368.
Wilton D R,Glisson AW,Schaubert D H,et al.Potentials integrals foruniform and linear source distributions on polygonal and polyhedral domains[J].IEEE Transactions on Antennas and Propagation,1984,32(3):276-281.

Claims (7)

1.一种重力正演加速方法，其特征在于，包括以下步骤：
步骤1、区域离散与八叉树建立：
首先，确定三维地质体Ω中的点的坐标和密度，并设计观测点的个数与坐标；三维地质体Ω中的点称为源点；
将三维地质体Ω剖分成多个四面体单元Ti,i＝1,2,...,N，其中N为四面体单元的个数；那么，有
然后，分别对观测点和三维地质体Ω建立对应的八叉树结构，分别称为：观测点八叉树和源八叉树；
对观测点建立对应的八叉树结构：构建一个包含所有观测点在内的六面体单元，称为最高层细胞，即第0层细胞，然后以第0层细胞为母细胞，将其分成8个子细胞，记为第1层细胞；继续按上述方式对每个子细胞进行剖分，后续每剖分一层得到的细胞依次记为第i层细胞，其中i＝2,3,4…,p，直至每个细胞包含的观测点数目少于设定的个数为止，并删除所有不含有观测点的细胞；最后一层的细胞称为叶子细胞；
同样对三维地质体Ω建立对应的八叉树结构：构建一个包含三维地质体Ω在内的六面体单元，称为最高层细胞，即第0层细胞，然后以第0层细胞为母细胞，将其分成8个子细胞，记为第1层细胞；继续按上述方式对每个子细胞进行剖分，后续每剖分一层得到的子细胞依次记为第j层细胞，其中j＝2,3,4…,q，直到每个细胞里包含的四面体单元数少于设定的个数为止，并删除所有不含有四面体单元的细胞；最后一层的子细胞称为叶子细胞；
步骤2、将三维地质体Ω划分为近区与远区：
对于每一个观测点r，若三维地质体Ω中某个四面体单元的中心点到观测点r的距离小于某一设定的数值ε，则认为该四面体单元所在区域为近区，记为Tnear；其他区域为远区，记为Tfar；有Ω＝Tnear+Tfar；
步骤3、分别计算近区与远区重力参数；
近区重力参数包括在观测点r处产生的近区重力位、近区重力梯度和近区重力梯度张量，分别记为和远区重力参数包括在观测点r处产生的远区重力位、远区重力梯度和远区重力梯度张量，分别记为和
首先，根据步骤2中划分的区域，将以下重力参数计算公式中三维地质体Ω的积分分解为近区的积分和远区的积分；
<mrow> <mi>&amp;phi;</mi> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>G</mi> <munder> <mrow> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> </mrow> <mi>&amp;Omega;</mi> </munder> <mfrac> <mrow> <mi>&amp;rho;</mi> <mrow> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> </mrow> </mrow> <mi>R</mi> </mfrac> <mi>d</mi> <mi>v</mi> </mrow>
<mrow> <mi>g</mi> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>G</mi> <munder> <mrow> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> </mrow> <mi>&amp;Omega;</mi> </munder> <mo>&amp;dtri;</mo> <mfrac> <mrow> <mi>&amp;rho;</mi> <mrow> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> </mrow> </mrow> <mi>R</mi> </mfrac> <mi>d</mi> <mi>v</mi> </mrow>
<mrow> <msub> <mi>T</mi> <mi>g</mi> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>G</mi> <munder> <mrow> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> </mrow> <mi>&amp;Omega;</mi> </munder> <mo>&amp;dtri;</mo> <mo>&amp;dtri;</mo> <mfrac> <mrow> <mi>&amp;rho;</mi> <mrow> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> </mrow> </mrow> <mi>R</mi> </mfrac> <mi>d</mi> <mi>v</mi> </mrow>
其中，为梯度运算符，R＝|r-r′|为观测点r到源点r′的距离；dv表示体积分ρ(r′)为三维地质体Ω在源点r′处的密度，G为引力常数，G＝6.674×10-11m3kg-1s-2；
然后，基于近区所包含的四面体单元的重力积分无奇异性解析表达式直接计算近区的积分，得到近区重力参数；并基于观测点八叉树和源八叉树，采用自适应快速多极展开算法计算远区的积分，得到远区重力参数；
步骤4、将步骤3中计算得到的近区重力参数与远区重力参数求和，得到观测点重力参数，包括在观测点r处产生的重力位φ(r)、重力梯度g(r)和重力梯度张量Tg(r)，其计算公式分别为：
<mrow> <mtable> <mtr> <mtd> <mrow> <mi>&amp;phi;</mi> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>&amp;phi;</mi> <msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>n</mi> <mi>e</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </msub> <mo>+</mo> <mi>&amp;phi;</mi> <msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>f</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>g</mi> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>g</mi> <msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>n</mi> <mi>e</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </msub> <mo>+</mo> <mi>g</mi> <msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>f</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>T</mi> <mi>g</mi> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>T</mi> <mi>g</mi> </msub> <msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>n</mi> <mi>e</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </msub> <mo>+</mo> <msub> <mi>T</mi> <mi>g</mi> </msub> <msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>f</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </msub> </mrow> </mtd> </mtr> </mtable> <mo>.</mo> </mrow>
2.根据权利要求1所述重力正演加速方法，其特征在于：设rc p为源点r′所在叶子细胞的中心点，满足条件|r-rc p|>|r′-rc p|；
首先，直接通过以下多极距公式计算rc p处的多极距Mnm(φ,r′,rc p)：
<mrow> <msub> <mi>M</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>,</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msup> <msub> <mi>r</mi> <mi>c</mi> </msub> <mi>p</mi> </msup> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mrow> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>f</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </munder> <msub> <mi>R</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>-</mo> <msup> <msub> <mi>r</mi> <mi>c</mi> </msub> <mi>p</mi> </msup> <mo>)</mo> </mrow> <mi>&amp;rho;</mi> <mrow> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> </mrow> <mi>d</mi> <mi>v</mi> <mo>;</mo> </mrow>
然后，通过以下公式计算远区重力位：
<mrow> <mi>&amp;phi;</mi> <msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>f</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </msub> <mo>=</mo> <mi>G</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>0</mn> </mrow> <mi>&amp;infin;</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mo>-</mo> <mi>n</mi> </mrow> <mi>n</mi> </munderover> <mover> <msub> <mi>S</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mo>&amp;OverBar;</mo> </mover> <mrow> <mo>(</mo> <mi>r</mi> <mo>-</mo> <msup> <msub> <mi>r</mi> <mi>c</mi> </msub> <mi>p</mi> </msup> <mo>)</mo> </mrow> <msub> <mi>M</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>,</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msup> <msub> <mi>r</mi> <mi>c</mi> </msub> <mi>p</mi> </msup> <mo>)</mo> </mrow> <mo>,</mo> <mo>|</mo> <mrow> <mi>r</mi> <mo>-</mo> <msup> <msub> <mi>r</mi> <mi>c</mi> </msub> <mi>p</mi> </msup> </mrow> <mo>|</mo> <mo>&gt;</mo> <mo>|</mo> <mrow> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>-</mo> <msup> <msub> <mi>r</mi> <mi>c</mi> </msub> <mi>p</mi> </msup> </mrow> <mo>|</mo> </mrow>
其中，Rnm和Snm分别为已知的常规的和奇异的球谐函数；为Snm的复共轭。
3.根据权利要求2所述重力正演加速方法，其特征在于：设点rc i为源八叉树第i层细胞的中心点，点rc i+1为源八叉树第i+1层细胞的中心点，点rc i所在细胞为点rc i+1所在细胞的母细胞；其中i＝0,1,2…,p-1；
首先，通过点rc i+1处的多极距与以下M2M转换公式计算点rc i处的多极距
<mrow> <msub> <mi>M</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>,</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msup> <msub> <mi>r</mi> <mi>c</mi> </msub> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <msup> <mi>n</mi> <mo>&amp;prime;</mo> </msup> <mo>=</mo> <mn>0</mn> </mrow> <mi>n</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <msup> <mi>m</mi> <mo>&amp;prime;</mo> </msup> <mo>=</mo> <mo>-</mo> <msup> <mi>n</mi> <mo>&amp;prime;</mo> </msup> </mrow> <msup> <mi>n</mi> <mo>&amp;prime;</mo> </msup> </munderover> <msub> <mi>R</mi> <mrow> <msup> <mi>n</mi> <mo>&amp;prime;</mo> </msup> <msup> <mi>m</mi> <mo>&amp;prime;</mo> </msup> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>-</mo> <msup> <msub> <mi>r</mi> <mi>c</mi> </msub> <mi>i</mi> </msup> <mo>)</mo> </mrow> <msub> <mi>M</mi> <mrow> <mi>n</mi> <mo>-</mo> <msup> <mi>n</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <mi>m</mi> <mo>-</mo> <msup> <mi>m</mi> <mo>&amp;prime;</mo> </msup> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>,</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>i</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>)</mo> </mrow> </mrow>
然后，通过以下公式计算远区重力位：
<mrow> <mi>&amp;phi;</mi> <msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>f</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </msub> <mo>=</mo> <mi>G</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>0</mn> </mrow> <mi>&amp;infin;</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mo>-</mo> <mi>n</mi> </mrow> <mi>n</mi> </munderover> <mover> <msub> <mi>S</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mo>&amp;OverBar;</mo> </mover> <mrow> <mo>(</mo> <mi>r</mi> <mo>-</mo> <msup> <msub> <mi>r</mi> <mi>c</mi> </msub> <mi>i</mi> </msup> <mo>)</mo> </mrow> <msub> <mi>M</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>,</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msup> <msub> <mi>r</mi> <mi>c</mi> </msub> <mi>i</mi> </msup> <mo>)</mo> </mrow> <mo>.</mo> </mrow>
4.根据权利要求3所述重力正演加速方法，其特征在于：设为观测点八叉树中最高层细胞的中心点；
首先，直接通过以下积分公式或M2L转换公式计算得到处的局部系数
积分公式为：
<mrow> <msub> <mi>L</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>,</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mrow> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>f</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </munder> <mover> <msub> <mi>S</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mo>&amp;OverBar;</mo> </mover> <mrow> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>-</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> <mi>&amp;rho;</mi> <mrow> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> </mrow> <mi>d</mi> <mi>v</mi> <mo>;</mo> </mrow>
M2L转换公式为：
<mrow> <msub> <mi>L</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>,</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <msup> <mi>n</mi> <mo>&amp;prime;</mo> </msup> <mo>=</mo> <mn>0</mn> </mrow> <mi>&amp;infin;</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <msup> <mi>m</mi> <mo>&amp;prime;</mo> </msup> <mo>=</mo> <mo>-</mo> <msup> <mi>n</mi> <mo>&amp;prime;</mo> </msup> </mrow> <msup> <mi>n</mi> <mo>&amp;prime;</mo> </msup> </munderover> <msup> <mrow> <mo>(</mo> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mi>n</mi> </msup> <mover> <msub> <mi>S</mi> <mrow> <mi>n</mi> <mo>+</mo> <msup> <mi>n</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <mi>m</mi> <mo>+</mo> <msup> <mi>m</mi> <mo>&amp;prime;</mo> </msup> </mrow> </msub> <mo>&amp;OverBar;</mo> </mover> <mrow> <mo>(</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mn>0</mn> </msubsup> <mo>)</mo> </mrow> <msub> <mi>M</mi> <mrow> <msup> <mi>n</mi> <mo>&amp;prime;</mo> </msup> <msup> <mi>m</mi> <mo>&amp;prime;</mo> </msup> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>,</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mn>0</mn> </msubsup> <mo>)</mo> </mrow> <mo>,</mo> <mo>|</mo> <mrow> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mn>0</mn> </msubsup> </mrow> <mo>|</mo> <mo>&gt;</mo> <mo>|</mo> <mrow> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>-</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mn>0</mn> </msubsup> </mrow> <mo>|</mo> </mrow>
其中，为处的多极距，通过M2M转换公式递推得到；
然后，通过以下公式计算远区重力位：
<mrow> <mi>&amp;phi;</mi> <msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>f</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </msub> <mo>=</mo> <mi>G</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>0</mn> </mrow> <mi>&amp;infin;</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mo>-</mo> <mi>n</mi> </mrow> <mi>n</mi> </munderover> <msub> <mi>R</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>-</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> <msub> <mi>L</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>,</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>,</mo> <mo>|</mo> <mrow> <mi>r</mi> <mo>-</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> </mrow> </msubsup> </mrow> <mo>|</mo> <mo>&lt;</mo> <mo>|</mo> <mrow> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>-</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> </mrow> </msubsup> </mrow> <mo>|</mo> <mo>.</mo> </mrow>
5.根据权利要求4所述重力正演加速方法，其特征在于：设点为观测点八叉树第j+1层细胞的中心点，点为观测点八叉树第j层细胞的中心点，点所在细胞为点所在细胞的母细胞；j＝0,1,2…,q-1；
首先，通过点处的局部系数和以下L2L转换公式计算点处的局部系数
<mrow> <msub> <mi>L</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>,</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <msup> <mi>n</mi> <mo>&amp;prime;</mo> </msup> <mo>=</mo> <mi>n</mi> </mrow> <mi>&amp;infin;</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <msup> <mi>m</mi> <mo>&amp;prime;</mo> </msup> <mo>=</mo> <mo>-</mo> <msup> <mi>n</mi> <mo>&amp;prime;</mo> </msup> </mrow> <msup> <mi>n</mi> <mo>&amp;prime;</mo> </msup> </munderover> <msub> <mi>R</mi> <mrow> <msup> <mi>n</mi> <mo>&amp;prime;</mo> </msup> <mo>-</mo> <mi>n</mi> <mo>,</mo> <msup> <mi>m</mi> <mo>&amp;prime;</mo> </msup> <mo>-</mo> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> <msub> <mi>L</mi> <mrow> <msup> <mi>n</mi> <mo>&amp;prime;</mo> </msup> <msup> <mi>m</mi> <mo>&amp;prime;</mo> </msup> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>,</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> </mrow>
然后，通过以下公式计算远区重力位：
<mrow> <mi>&amp;phi;</mi> <msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>f</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </msub> <mo>=</mo> <mi>G</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>0</mn> </mrow> <mi>&amp;infin;</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mo>-</mo> <mi>n</mi> </mrow> <mi>n</mi> </munderover> <msub> <mi>R</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>-</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> <msub> <mi>L</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>,</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>.</mo> </mrow>
6.根据权利要求5所述重力正演加速方法，其特征在于：通过以下公式计算远区重力场和远区重力梯度张量：
<mrow> <mtable> <mtr> <mtd> <mrow> <mi>g</mi> <msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>f</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </msub> <mo>=</mo> <mi>G</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>0</mn> </mrow> <mi>&amp;infin;</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mo>-</mo> <mi>n</mi> </mrow> <mi>n</mi> </munderover> <mo>&amp;dtri;</mo> <msub> <mi>R</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>-</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> <msub> <mi>L</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>,</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>T</mi> <mi>g</mi> </msub> <msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>f</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </msub> <mo>=</mo> <mi>G</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>0</mn> </mrow> <mi>&amp;infin;</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mo>-</mo> <mi>n</mi> </mrow> <mi>n</mi> </munderover> <mo>&amp;dtri;</mo> <mo>&amp;dtri;</mo> <msub> <mi>R</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>-</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> <msub> <mi>L</mi> <mrow> <mi>n</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;phi;</mi> <mo>,</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msubsup> <mi>r</mi> <mi>c</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> <mo>.</mo> </mrow>
7.根据权利要求1～6中任一项所述的重力正演加速方法，其特征在于：采用基于近区所包含的四面体单元的重力积分无奇异性解析表达式直接计算近区重力参数的计算公式为：
<mrow> <mi>&amp;phi;</mi> <msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>n</mi> <mi>e</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </msub> <mo>=</mo> <mi>G</mi> <munder> <mrow> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>n</mi> <mi>e</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </munder> <mfrac> <mrow> <mi>&amp;rho;</mi> <mrow> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> </mrow> </mrow> <mi>R</mi> </mfrac> <mi>d</mi> <mi>v</mi> </mrow>
<mrow> <mi>g</mi> <msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>n</mi> <mi>e</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </msub> <mo>=</mo> <mi>G</mi> <munder> <mrow> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>n</mi> <mi>e</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </munder> <mo>&amp;dtri;</mo> <mfrac> <mrow> <mi>&amp;rho;</mi> <mrow> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> </mrow> </mrow> <mi>R</mi> </mfrac> <mi>d</mi> <mi>v</mi> </mrow>
<mrow> <msub> <mi>T</mi> <mi>g</mi> </msub> <msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>n</mi> <mi>e</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </msub> <mo>=</mo> <mi>G</mi> <munder> <mrow> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> <mo>&amp;Integral;</mo> </mrow> <msub> <mi>T</mi> <mrow> <mi>n</mi> <mi>e</mi> <mi>a</mi> <mi>r</mi> </mrow> </msub> </munder> <mo>&amp;dtri;</mo> <mo>&amp;dtri;</mo> <mfrac> <mrow> <mi>&amp;rho;</mi> <mrow> <mo>(</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> </mrow> </mrow> <mi>R</mi> </mfrac> <mi>d</mi> <mi>v</mi> </mrow>
其中，R＝1/|r-r′|。
CN201611252162.XA 2016-12-29 2016-12-29 一种重力正演加速方法 Active CN106646645B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201611252162.XA CN106646645B (zh) 2016-12-29 2016-12-29 一种重力正演加速方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201611252162.XA CN106646645B (zh) 2016-12-29 2016-12-29 一种重力正演加速方法

Publications (2)

Publication Number Publication Date
CN106646645A CN106646645A (zh) 2017-05-10
CN106646645B true CN106646645B (zh) 2018-01-19

Family

ID=58836767

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201611252162.XA Active CN106646645B (zh) 2016-12-29 2016-12-29 一种重力正演加速方法

Country Status (1)

Country Link
CN (1) CN106646645B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107561592B (zh) * 2017-09-11 2018-07-17 中南大学 一种密度为多项式的任意多面体的重力场正演方法
CN108984939A (zh) * 2018-07-30 2018-12-11 中南大学 基于3D Gauss-FFT的三维重力场正演方法
CN109444973B (zh) * 2018-11-06 2019-12-13 中南大学 一种球坐标系下重力正演加速方法
CN109375280B (zh) * 2018-12-10 2020-03-31 中南大学 一种球坐标系下重力场快速高精度正演方法
CN112965125B (zh) * 2021-02-08 2022-08-05 中国人民解放军92859部队 一种基于重力异常计算外部扰动重力东向分量的方法
CN113204057B (zh) * 2021-05-07 2022-07-12 湖南科技大学 一种基于多层级方法重磁快速反演方法
CN113238284B (zh) * 2021-05-07 2022-09-27 湖南科技大学 一种重磁快速正演方法
CN113109883B (zh) * 2021-05-26 2023-01-03 湖南科技大学 基于等参变换全球离散网格球坐标下卫星重力场正演方法
CN113642189B (zh) * 2021-08-25 2023-09-19 中南大学 一种基于积分解的重力梯度张量快速数值模拟方法和装置

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105334542A (zh) * 2015-10-23 2016-02-17 中南大学 任意密度分布复杂地质体重力场快速、高精度正演方法
CN105549106A (zh) * 2016-01-07 2016-05-04 中国科学院地质与地球物理研究所 一种重力多界面反演方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2014186479A2 (en) * 2013-05-15 2014-11-20 Schlumberger Canada Limited Geobody surface reconstruction

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105334542A (zh) * 2015-10-23 2016-02-17 中南大学 任意密度分布复杂地质体重力场快速、高精度正演方法
CN105549106A (zh) * 2016-01-07 2016-05-04 中国科学院地质与地球物理研究所 一种重力多界面反演方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
三维重力梯度快速正反演研究;陈涛;《中国优秀硕士学位论文全文数据库 基础科学辑》;20140815;全文 *
基于GPU并行的重力、重力梯度三维正演快速计算及反演策略;陈召曦;《地球物理学报》;20121231;第55卷(第12期);全文 *
基于OpenMP的重力张量并行正演;余涛;《物探化探计算技术》;20130731;第35卷(第4期);全文 *

Also Published As

Publication number Publication date
CN106646645A (zh) 2017-05-10

Similar Documents

Publication Publication Date Title
CN106646645B (zh) 一种重力正演加速方法
CN109375280B (zh) 一种球坐标系下重力场快速高精度正演方法
Praveen et al. Low cost PSO using metamodels and inexact pre-evaluation: Application to aerodynamic shape design
Han et al. Accuracy of the adaptive GRP scheme and the simulation of 2-D Riemann problems for compressible Euler equations
CN107561592B (zh) 一种密度为多项式的任意多面体的重力场正演方法
CN116738891B (zh) 一种增强飞行器流场模拟稳定性的lu-sgs改进方法
CN115292973B (zh) 一种任意采样的空间波数域三维磁场数值模拟方法及系统
Xia et al. OpenACC acceleration of an unstructured CFD solver based on a reconstructed discontinuous Galerkin method for compressible flows
CN106855635A (zh) 一种确定正交观测系统的面元属性信息的方法和装置
CN114722731B (zh) 一种拖体粘滞拖曳阻力预测方法及系统
CN116774292A (zh) 一种地震波走时确定方法、系统、电子设备及存储介质
Liang et al. Uncertainty quantification of geologic model parameters in 3D gravity inversion by Hessian-informed Markov chain Monte Carlo
CN106662665B (zh) 用于更快速的交错网格处理的重新排序的插值和卷积
CN111158059B (zh) 基于三次b样条函数的重力反演方法
CN110967778B (zh) 一种动态坐标系多面体剖分重力布格校正方法
CN115640720B (zh) 一种基于距离控制网格加密的自引力仿真方法
Mouallem et al. Implementation of the Novel Duo‐Grid in GFDL's FV3 Dynamical Core
CN111046542B (zh) 一种评估二十面体的le网格计算特性方法及离散方法
CN113267830A (zh) 基于非结构网格下二维重力梯度与地震数据联合反演方法
CN111950108B (zh) 二度变密度体重力梯度张量计算方法
CN108205161B (zh) 一种基于大数据计算的航空重力数据物性层析处理方法
CN112748471A (zh) 一种非结构化等效源的重磁数据延拓与转换方法
CN110568497A (zh) 一种复杂介质条件下地震初至波旅行时的精确求解方法
RU2611892C1 (ru) Способ трехмерного моделирования заданного гидрогеологического объекта, реализуемый в вычислительной системе
Qin et al. A method for calculating two-dimensional spatially extension distances and its clustering algorithm

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